Bifurcation and optimal control analysis of HIV/AIDS and COVID-19 co-infection model with numerical simulation

HIV/AIDS and COVID-19 co-infection is a common global health and socio-economic problem. In this paper, a mathematical model for the transmission dynamics of HIV/AIDS and COVID-19 co-infection that incorporates protection and treatment for the infected (and infectious) groups is formulated and analyzed. Firstly, we proved the non-negativity and boundedness of the co-infection model solutions, analyzed the single infection models steady states, calculated the basic reproduction numbers using next generation matrix approach and then investigated the existence and local stabilities of equilibriums using Routh-Hurwiz stability criteria. Then using the Center Manifold criteria to investigate the proposed model exhibited the phenomenon of backward bifurcation whenever its effective reproduction number is less than unity. Secondly, we incorporate time dependent optimal control strategies, using Pontryagin’s Maximum Principle to derive necessary conditions for the optimal control of the disease. Finally, we carried out numerical simulations for both the deterministic model and the model incorporating optimal controls and we found the results that the model solutions are converging to the model endemic equilibrium point whenever the model effective reproduction number is greater than unity, and also from numerical simulations of the optimal control problem applying the combinations of all the possible protection and treatment strategies together is the most effective strategy to drastically minimizing the transmission of the HIV/AIDS and COVID-19 co-infection in the community under consideration of the study.


Introduction
Infectious diseases are diagnostically proven illnesses caused by tiny microorganisms such as viruses, bacteria, fungi, and parasites and have been the leading causes of death throughout the world, for example; viruses cause both COVID-19 and HIV/AIDS infections [1][2][3].
Human immunodeficiency virus (HIV) is one of the most dangerous viruses that is spreading around the world. AIDS, or acquired immunodeficiency syndrome, is one of the most devastating epidemics in history, caused by HIV, which has been a worldwide epidemic since 1981 [4][5][6][7][8][9][10]. It remains a significant world health issue that impacts almost seventy million people worldwide and has been a significant cause of morbidity and mortality [11,12]. HIV is transmissible through sexual contact, needle sharing, and direct contact with virus-infected blood or other body fluids, as well as from mother to child during giving birth [10,[13][14][15]. In early December 2019, a coronavirus called COVID-19 was reported in Wuhan, China, with symptoms similar to pneumonia. According to reports, it is one of the most devastating infectious diseases caused by the novel coronavirus SARS-CoV-2, which has been a significant impact on the health, social, and economic integration of communities worldwide [16][17][18][19][20][21][22][23][24][25][26][27][28][29]. On March 11, 2020, the World Health Organization (WHO) confirmed it as a global pandemic, and on July 25, 2020, the world total number of COVID-19 infected individuals was 15,762,007, with 640,276 deaths [25,28,29]. It was suspected to be pneumonia or a common cold-like illness, with symptoms such as fatigue, alter in taste, fever, muscular pains, shortness of breath, ironical cough, and sore throat [25,27,30]. Despite massive efforts to reduce the virus's transmission and survivability, the death rate from COVID-19 remains high [15]. COVID-19 can be transmitted through sneezing or coughing droplets expelled from the human lungs, as well as when humans come into contact with contaminated dispatched materials [17,26,31]. Among the unfortunate aspects of the COVID-19 pandemic is that patients over the age of 60 are more likely to be infected than anyone below the age of 60 [31]. It is an extremely infectious contagious agent that has spread throughout most of the world's nations and has a significant impact on the global economy and public health [24,32]. COVID-19 infection may be more common in people with compromised immunity from other infections such as tuberculosis, HIV, pneumonia, and cholera [1,25,[33][34][35][36][37]. WHO unanimously implemented vaccination, quarantine, wearing face masks, hand washing with alcohol, and significant discrepancies as possible prevention and control strategies [26,27,31]. Symptomless and pre-symptomatic transmission, a low incidence or lack of dominant systemic symptoms such as fever, airborne transmission that may require a high infectious dose and super-spread events are the essential aspects of COVID-19 spreading that make it challenging to handle [16].
Mathematical modelling approaches have been crucial to provide basic frameworks in order to understand the transmission dynamics of infectious diseases [37]. Many scholars throughout the world have been formulated and analyzed mathematical models to investigate the transmission dynamics of different infectious diseases using ordinary differential equations approach like [2,9,15,17,19,22,23,[26][27][28][29]31,32,[45][46][47] using stochastic approach like [48], and using fractional order derivative approach like [1,5,49,50]. In the structure of this study, we have reviewed research papers that have been done on the transmission dynamics of different infectious diseases especially co-infections of HIV/AIDS and other infectious diseases. Teklu and Rao [14] constructed and examined HIV/AIDS and pneumonia co-infection model with control measures such as pneumonia vaccination and treatments of pneumonia and HIV/ AIDS infections. Hezam et al. [40], formulated a mathematical model for cholera and COVID-19 co-infection which describes the transmission dynamics of COVID-19 and cholera in Yemen. The model analysis examined four controlling measures such as social distancing, lockdown, the number of test kits to control the COVID-19 outbreak, and the number of susceptible individuals who can get CWTs for water purification. Anwar et al. [15], constructed a mathematical model on COVID-19 with the isolation controlling measure on the COVID-19 infected individuals throughout the community. Ahmed et al. [1] formulated and analyzed HIV and COVID-19 co-infection model with ABC-fractional operator approach to investigate an epidemic prediction of a combined HIV-COVID-19 co-infection model. Numerical simulations were carried out to justify that the disease will stabilize at a later stage when enough protection strategies are taken. Teklu and Terefe [3] analyze COVID-19 and syphilis co-dynamics model to investigate the impacts of intervention measures on the disease transmission.
Similarly, various Scholars have formulated and analyzed mathematical models with optimal control strategies to investigate the effect of prevention and control measures on HIV/ AIDS, COVID-19, HIV/AIDS and COVID-19 co-infection and other various infectious diseases transmission throughout nations in the world. For instance, Tchoumi et al. [37] proposed and investigated the co-dynamics of malaria and COVID-19 co-dynamics: with optimal control strategies. The numerical simulation results verifies the theoretical optimal control analysis and illustrates that using malaria and COVID-19 protection measures concurrently can help mitigate there transmission compared with applying single infections protection measures. Omame et al. [25] investigated a mathematical model for the dynamics of COVID-19 infection in order to assess the impacts of prior comorbidity on COVID-19 complications and COVID-19 reinfection with optimal control strategies. The authors recommended that the strategy that prevents COVID-19 infection by comorbid susceptible is the best cost-effective of all the other control strategies for the prevention of COVID-19. Ringa et al. [43] formulated and analyzed a mathematical model on HIV and COVID-19 co-infection with optimal control strategies. Their analysis suggested that COVID19 only prevention strategy is the most effective strategy and it averted about 10,500 new co-infection cases. Keno et al. [51] investigated an optimal control and cost effectiveness analysis of SIRS malaria disease model with temperature variability facto. Their result suggested that the combination of treatment of infected humans and insecticide spraying was proved to be the best efficient and least costly strategy to eradicate the disease. Keno et al. [52] investigated a mathematical model with optimal control strategies for malaria transmission with role of climate variability. Their result suggested that the combination of treated bed net and treatment is the most optimal and least-cost strategy to minimize the malaria. Goudiaby et al. [39] formulated and analyzed a COVID-19 and tuberculosis co-dynamics model with optimal control strategies. They suggested that COVID-19 prevention, treatment and control of co-infection yields a better outcome in terms of the number of COVID-19 cases prevented at a lower percentage of the total cost of this strategy. Asamoah et al. [53] constructed a mathematical model on COVID-19 to investigate optimal control strategies and comprehensive cost-effectiveness. Okosun et al. [54] formulated a mathematical model on HIV/AIDS to investigate the impact of optimal control on the treatment of HIV/ AIDS and screening of unaware invectives. Their analysis recommended that the combination of all the control strategies is the most cost-effective strategy. Furthermore, notice that optimal control modeling and cost-effectiveness analysis model have been applied in recent infectious diseases models like [55,56].
As we observed from review of literatures done by various epidemiology and medical scholars, HIV/AIDS and COVID-19 co-infection is a public health concern especially in developing nations of the world. The main purpose of this paper is to investigate the impacts of COVID-19 protection with quarantine, COVID-19 treatment, HIV protection and HIV treatment prevention and controlling strategies on the transmission dynamics of HIV/AIDS and COVID-19 co-infection in the community with mathematical modelling approach. We have reviewed literatures [1,43] invested much effort in studying HIV/AIDS and COVID-19 co-infection, but did not considered COVID-19 protection with quarantine, COVID-19 treatment, HIV/AIDS protection, and HIV/AIDS treatment as prevention and control strategies simultaneously in a single model formulation which motivates us to undertake this study and fill the gap.

Basic frameworks of the model
In this paper, we partitioned the total human population at a given time t denoted by N(t), into eleven mutually-exclusive classes depending on their infection status: susceptible class to both COVID-19 and HIV S(t)), COVID-19 protection by quarantine class (C q (t)), HIV protected (such as by using condom, limit sexual partners, creating awareness etc.) class (H p (t)), COVID-19 protection by vaccination class (C v (t)), COVID-19 mono-infection class (C i (t)), HIV unaware mono-infection class (H u (t)), HIV aware mono-infection class (H a (t)), HIV unaware and COVID-19 co-infection class (M u (t)), HIV aware and COVID-19 co-infection class (M a (t)), COVID-19 recovery class (R(t)), and HIV aware treatment class (H t (t)) so that; Since HIV is a chronic infectious disease the susceptible individuals acquires HIV infection at the standard incidence rate given by where ρ 3 � ρ 2 � ρ 1 � 1 are the modification parameters that increase infectivity and β 1 is the HIV transmission rate. Since COVID-19 is a very acute infection the susceptible individuals acquires COVID-19 infection at the mass action incidence rate as stated in [50,51,54].
where ω 2 > ω 1 > 1 are the modification parameters that increase infectivity and β 2 is the COVID-19 transmission rate. Additional model assumptions • k 1 , k 2 , k 3 , and k 4 where k 4 = 1 − k 1 − k 2 − k 3 are portions of the number of recruited individuals those are entering to the susceptible class, the COVID-19 protected class, the HIV protected class and the COVID-19 vaccination class respectively.
• The susceptible class is increased by individuals from the COVID-19 vaccinated class in which those individuals who are vaccinated against COVID-19 but did not respond to vaccination with waning rate of ρ and from COVID-19 recovery with treatment class who develop their temporary immunity by the rate η.
• COVID-19 vaccine is may not be 100% efficient, so vaccinated individuals also have a chance of being infected with portion ε of the serotype not covered by the vaccine where 0 � ε < 1.
• 0 < υ � 1 is the modification parameter such that COVID-19 infected individual is less susceptible to HIV infection than a susceptible individuals due to morbidity.
• There is screening and testing mechanisms for the previous and current status in each class.
• The human population distribution is homogeneous in each class.
• HIV treated individuals do not transmit infection to others due to awareness.
• Population of human being is variable.

PLOS ONE
Bifurcation analysis and optimal control for COVID-19 and HIV/AIDS co-infection • There is no dual-infection transmission simultaneously.
• No vertical HIV transmission.
In this section using parameters given in Table 1, model variables given in Table 2, and the model basic frame work, and assumptions given in (2.1), the schematic diagram for the transmission dynamics of HIV/AIDS and COVID-19 co-infection is given by   with the corresponding initial conditions The sum of all the differential equations in (3) is

The basic qualitative properties of the model (3)
The COVID-19 and HIV/AIDS co-infection model given in Eq (3) Theorem 1 (Positivity of the model solutions) Let us given the initial data in Eq (4) then the solutions Proof: Let us consider We have to show that and H t (t) > 0}. Now since the entire co-infection model state variables are positive and all the state variables are continuous, we can justify that τ > 0. If τ = +1, then non-negativity holds. But, if 0 < τ < +1 we will have S(τ) = 0 or C q (τ) = 0 or H p (τ) = 0 or C v (τ) = 0 or C i (τ) = 0 or H u (τ) = or H a (τ) = 0 or M u (τ) = 0 or M a (τ) = 0 or R(τ) = 0 or H t (τ) = 0.
Here from the first equation of the COVID-19 and HIV/AIDS co-infection model (3) we have got and integrate using method of integrating factor we have determined the constant value and from the meaning of τ, the solutions C q (t) > 0, H p (t) > 0, C v (t) > 0, R(t) > 0. Moreover, the exponential function is always positive, then the solution S(τ) > 0 hence S(τ) 6 ¼ 0. Thus following the same procedure for τ = +1, all the solutions of the COVID-19 and HIV/AIDS coinfection system (3) are non-negative.

Theorem 2 (The invariant region):
All the feasible positive solutions of the co-infection model (3) are bounded in the region (6).
þ is an arbitrary non-negative solution of the system (3) with initial conditions given in Eq (4). Now adding all the differential equations given in Eq (3) we have got the derivative of the total population N which is given in Eq (5) as Then by ignoring the infections we have determined that _ N � D À mN and using separation of variables whenever t ! 1, we have obtained that 0 � N � D m . Hence, all the positive feasible solutions of the co-infection model (3) entering in to the region given in Eq (6).
Note: Since the model (3) solutions are both positive and bounded in the region (6) the HIV/AIDS and COVID-19 co-infection model (3) is both mathematically and biologically meaning full [45,47,57], then we can consider the two mono-infection models, namely; HIV mono-infection and COVID-19 mono-infection models. This is fundamental for the analysis of the COVID-19 and HIV/AIDS co-infection model.

Analytical result of the models
Before analyzing the HIV/AIDS and COVID-19 co-infection model given in Eq (3), it is very crucial to gain some basic backgrounds about the COVID-19 and HIV/AIDS mono-infection models.

Mathematical analysis of HIV/AIDS mono-infection model
In this subsection we assume there is no COVID-19 infection in the community i.e. C q =, C q = C i = M u = M a = R = 0 in (3) then the HIV/AIDS sub-model is given by , and the HIV sub-model force of infection given by In a similar manner of the full co-infection model (3) in it is sufficient to consider the dynamics of the sub-model (7) in O 1 as biologically and mathematically well-posed.

Disease-free equilibrium point of HIV mono-infection model (7) local stability
The disease-free equilibrium point of the HIV mono-infection system in (7) is obtained by making its right-hand side is equal to zero and setting the infected classes and treatment class to zero as H u = H a = H t = 0 which yields, Hence the disease-free equilibrium point is given by E 0 The local stability of the HIV mono-infection model (7) disease-free equilibrium point is examined by its effective reproduction number denoted by R HM , which is calculated by using the next generation operator method determined by Van den Driesch and Warmouth stated in [2]. Applying the method stated in [29], the transmission matrix F and the transition matrix V i.e., for the new infection and the remaining transfer respectively, are given by After some computations we have determined that and Then, the effective reproduction number of the HIV mono-infection model (7) is defined as the largest eigenvalue in magnitude of the next generation matrix, FV −1 given by The value R HM is defined as the total average number of secondary HIV unaware and HIV aware infection cases acquired from a typical HIV unaware or HIV aware individual during his/her effective infectious period in a susceptible population. The threshold result R HM is the effective reproduction number for HIV mono-infection.
Theorem 3: The disease-free equilibrium point of the HIV mono-infection model given in Eq (7) is locally asymptotically stable (LAS) if R HM < 1, and it is unstable if R HM > 1.
Proof: The local stability of the disease-free equilibrium point of HIV mono-infection model (7) is evaluated by applying the Routh-Hurwitz stability criteria stated in [52].
The Jacobian matrix of the HIV mono-infection model given in Eq (7) at the disease-free equilibrium point E 0 HM is given by : Then the corresponding characteristic equation of the Jacobian matrix JðE 0 HM Þ is given by On Eq (8) we applied Routh-Hurwitz stability criteria stated in [47] and we have determined that both eigenvalues are negative if R HM < 1. Furthermore, we can conclude that the disease-free equilibrium point of the model (7) is locally asymptotically stable whenever R HM < 1 since all the eigenvalues are negative when R HM < 1. The biological meaning of Theorem 3 can be stated as HIV infection can be eradicated from the population (whenever R HM < 1) if the initial size of the sub-populations of the HIV mono-infection model given in Eq (7) is in the basin of attraction of the disease-free equilibrium point E 0 HM .

Existence of HIV mono-infection endemic equilibrium point(s)
be an arbitrary endemic equilibrium point of the HIV monoinfection model (7) which can be determined by making the right hand side of Eq (7) as zero.
The after a number of steps of computations we have got where m 1 = (α 2 + μ), m 2 = (θ + μ + d 2 ), and m 3 = (γ + d 3 + μ). Now substitute H * u and H * a given in Eq (9) in to the HIV/AIDS force of infection Then we have the result Then the non-zero solution of (10) is l * H ¼ m 4 À m 5 m 6 . Therefore, the required non-zero solution (force of infection is obtained as l * H ¼

COVID-19 sub-model analysis
The corresponding COVID-19 sub-model of the system (3) is determined by making H p = H a = H u = M u = M a = H t = 0, and it is given by Here like the full model (3) and the HIV/AIDS sub-model (7) in the , it is sufficient to consider the dynamics of model (11) in O 2 be both biologically and mathematically meaningful.

Local stability of COVID-19 mono-infection model (11)
Disease-free equilibrium. Disease-free equilibrium point of the COVID-19 mono-infection model (11) is obtained by making its right-hand side as zero and setting the infected class and recovered with treatment class to zero as C i = R = 0 and after some simple steps of calculations we have determined that Hence the COVID-19 mono-infection model (11) disease-free equilibrium point is given by Here we are applying the Van Den Driesch and Warmouth next-generation matrix approach stated in [2] to determine the COVID-19 mono-infection model (11) effective reproduction number R C . After long computations, we have determined the transmission matrix given by and the transition matrix given by Then using Mathematica we have determined as and The characteristic equation of the matrix FV −1 is Then the spectral radius (effective reproduction number

Theorem 5:
The Disease-free equilibrium point E 0 CM of the COVID-19 mono-infection model (11) is locally asymptotically stable if R C < 1 otherwise unstable.

Proof:
The local stability of the disease-free equilibrium of the system (11) can be studied from its Jacobian matrix and Routh-Hurwitz stability criteria. The Jacobian matrix of the dynamical system at the disease-free equilibrium point is given by : Then the characteristic equation of the above Jacobian matrix is given by Therefore, since all the eigenvalues of the characteristics polynomials of the system (11) are negative if R C < 1 the disease-free equilibrium point of the COVID-19 mono-infection model (11) is locally asymptotically stable.

Existence of endemic equilibrium point (s) of the COVID-19 mono-infection model.
Before checking the global stability of the disease-free equilibrium point of the COVID-19 mono-infection model (11), we shall find the possible number of endemic equilibrium point(s) of the model (11). Let E * C ¼ ðS * ; C * q ; C * v ; C * i ; R * Þ be the endemic equilibrium point of COVID-19 mono-infection and l * C ¼ b 2 C * i be the COVID-19 mono-infection mass action incidence rate ("force of infection") at the equilibrium point. To find equilibrium point (s) for which COVID-19 mono-infection is endemic in the population, the equations are solved in terms of l * C ¼ b 2 C * i at an endemic equilibrium point. Now setting the right-hand sides of the equations of the model to zero (at steady state) gives and Then we have substituted in the COVID-19 force of infection given by l * C ¼ b 2 R * we have got the non-zero solution of l * C is obtained from the cubic equation where It can be seen from and (13) that c 3 > 0 (since the entire model parameters are nonnegative). Furthermore, c 0 > 0 whenever R C < 1. Thus, the number of possible positive real roots the polynomial (12) can have depends on the signs of c 1 , and c 2 . This can be analyzed using the Descartes' rule of signs on the cubic . Hence, the following results are established. (ii) c 1 < 0 and c 2 < 0.
(b). more than one endemic equilibrium point if R C > 1 either of the following holds.
Here, item (c) shows the happening of the backward bifurcation in the model (11) i.e., the locally asymptotically stable disease-free equilibrium point co-exists with a locally asymptotically stable endemic equilibrium point if R C < 1; examples of the existence of backward bifurcation phenomenon in mathematical epidemiological models, and the causes, can be seen in [8,17,26,31,[58][59][60]. The epidemiological consequence is that the classical epidemiological requirement of having the reproduction number R C to be less than one, even though necessary, is not sufficient for the effective control of the disease. The existence of the backward bifurcation phenomenon in sub-model (11) is now explored.
Theorem 7: The COVID-19 mono-infection model (11) exhibits backward bifurcation at R C ¼ 1 whenever the inequality D 2 > D 1 holds, where ÞðmþZÞ . In this section, we have used the center manifold theory stated in [60] to ascertain the local asymptotic stability of the endemic equilibrium due to the convolution of the first approach (eigenvalues of the Jacobian). To make use of the center manifold theory, the following change of variables is made by symbolizing S = x 1 , C p = x 2 , C v = x 3 , C i = x 4 and R = x 5 such that N 2 = x 1 + x 2 + x 3 + x 4 + x 5 . Furthermore, by using vector notation X = (x 1 , x 2 , x 3 , x 4 , x 5 ) T , the COVID-19 mono-infection model (11) can be written in the form dX with λ C = β 2 x 4 then the method entails evaluating the Jacobian of the system (14) Consider, R C ¼ 1 and suppose that β 2 = β* is chosen as a bifurcation parameter. From Solving for β 2 we have got b 2 ¼ b * ¼ mða 1 þmÞðrþmÞðmþd 1 þkÞÞ k 1 Dða 1 þmÞðrþmÞþa 1 k 2 DðrþmÞþk 4 Dða 1 þmÞðrþmεÞ .
After some steps of the calculation we have determined the eigenvalues of J β* as λ 1 = −μ, λ 2 = −(α 1 + μ) or or λ 3 = −(ρ + μ) or λ 4 = 0 or λ 5 = −(μ + η). It follows that the Jacobian J E 0 CM À � of Eq (14) at the disease-free equilibrium with β 2 = β*, denoted by J β* , has a simple zero eigenvalue with all the remaining eigenvalues have negative real part. Hence, Theorem 2 of Castillo-Chavez and Song stated in [60] can be used to analyze the dynamics of the model to show that the model (11) undergoes backward bifurcation at R C ¼ 1.
Eigenvectors of J β* : For the case R C ¼ 1, it can be shown that the Jacobian of the system (14) at β 2 = β* (denoted by J β* ) has a right eigenvectors associated with the zero eigenvalue given by u = (u 1 , u 2 , u 3 , u 4 , u 5 Then solving Eq (15) the right eigenvectors associated with the zero eigenvalue are given by Similarly, the left eigenvector associated with the zero eigenvalues at β 2 = β* given by ÞðmþZÞ . Thus, the bifurcation coefficient a is positive whenever

PLOS ONE
Hence, from the theory of Castillo-Chavez and Song stated in [60] the COVID-19 monoinfection model (11) exhibits a phenomenon of backward bifurcation at R C ¼ 1 and whenever D 2 > D 1 .
The diagram representation of this bifurcation is given in Fig 2 below. Fig 2 shows the appearance of backward bifurcation, which results in the coexistence of several equilibrium points. In such a case, the common conditions of disease eradication such as making R C < 1 will not work, and the initial number of infected persons also plays a crucial role.

Effective reproduction number of the co-infection model.
The effective reproduction number of the dynamical system (3) by applying the next generation operator method is the largest (dominant) eigenvalue (spectral radius) of the matrix: where F i is the rate of appearance of new infection in compartment i, v i is the transfer of infections from one compartment i to another, and E 0 is the disease-free equilibrium point. After ; where : Applying Mathematica we have determined as : After some computations and simplifications we have determined the dominant eigenvalue in magnitude of the matrix FV −1 which is the HIV/AIDS and COVID-19 co-infection effective reproduction number given by

Global asymptotic stability of disease-free equilibrium point.
In this sub-section we have used the method derived by Castillo-Chavez et al. and stated in reference [61] to look into the global asymptotic stability (GAS) of the co-infection model (3) disease-free equilibrium point. We mention two requirements that, if satisfied, also ensure the disease-free equilibrium is globally asymptotically stable. Then the new system (3) is rewritten as: where C ¼ ðS; C q ; H p ; C v Þ 2 R 4 denotes the number of uninfected components and z ¼ ðC i ; H u ; H a ; M u ; M a ; R; H t Þ 2 R 7 denotes the number of infected components. P 0 = (C 0 , 0), denotes the disease-free equilibrium point of the system. The following requirements must be satisfied to ensure the globally asymptotic stability: (H 1 ) For dC dt ¼ F C; 0 ð Þ, P 0 is globally asymptotically stable.
where C represents the number of non-infectious compartments and Y represents the number of infectious compartments. And

PLOS ONE
Bifurcation analysis and optimal control for COVID-19 and HIV/AIDS co-infection It is clear from the above discussion, that, b G C; Υ ð Þ≱0. Hence by the same reason given by results in reference [38], the disease-free equilibrium point may not be globally asymptotically stable.

Analysis of the optimal control strategy
In this section, we provide a thorough qualitative analysis of the time-dependent HIV/AIDS and COVID-19 co-infection model (3). The Pontryagin's Maximum Principle stated in literatures [25,43,51,52,55] is used to describe this analysis, with the aim of minimizing the HIV/ AIDS infection aware individuals denoted by H a , the COVID-19 infected individuals denoted by C i and the total HIV/AIDS and COVID-19 co-infected individuals denoted by M u + M a . In the case of time-dependent optimal control, we employ Pontryagin's Maximum Principle to derive the necessary conditions for diseases control mechanisms. After incorporating the controls into the HIV/AIDS and COVID-19 co-infection transmission model (3), the optimal control problem is as follows: with the corresponding initial conditions and 0 � u 1 t ð Þ � 1 represents HIV/AIDS infection protective control, 0 � u 2 t ð Þ � 1 represents the COVID-19 infections protective control using quarantine, 0 � u 3 t ð Þ � 1 represents the COVID-19 infection treatment control, and 0 � u 4 t ð Þ � 1 represents the HIV/AIDS treatment control.
The objective is to find the optimal control values u * ¼ u * 1 ; u * 2 ; u * 3 ; u * 4 À � of the controls u ¼ u 1 ; u 2 ; u 3 ; u 4 ð Þ such that the associated state trajectories are solution of the optimal control system (17) in the intervention time interval [0, T f ] with initial conditions as given in (18) and minimize the objective functional given by where the coefficients w 1 ; w 2 ; w 3 ; and w 4 are positive weight constants and B 1 2 ; B 2 2 ; B 3 2 and B 4 2 are the measure of relative costs of interventions associated with the controls u 1 ; u 2 ; u 3 and u 4 , respectively, and also balances the units of integrand. In the cost functional, the term w 1 A 5 refer to the cost related to COVID-19 infected class, the term w 2 H u refer to the cost related to individuals mono-infected with HIV and aware, the term w 3 A 8 refer to the cost related to co-infected individuals unaware of HIV infection and the term w 4 M a refer to the cost related to co-infected individuals aware of HIV infection.
, measures the current cost at time t. The set of admissible Lebesgue measurable control functions is defined by More precisely, we seek an optimal control pair Theorem 9 (Existence Theorem): There exists an optimal control u * ¼ u * 1 ; u * 2 ; u * 3 ; u * 4 À � in O u and a corresponding solution vector S * ; the optimal control dynamical system (17) with the initial values (18) We utilize Pontryagin's Maximal principle stated in literatures [51,52,55], to determine the prerequisites for the optimal control model (17). The optimal control problem (17) and (19) defined Hamiltonian (H) function is expressed as where G i stands for the i th state variable equation and λ 1 (t), λ 2 (t), λ 3 (t), λ 4 (t), λ 5 (t), λ 6 (t), λ 7 (t), λ 8 (t), λ 9 (t), λ 10 (t) and λ 11 (t) are adjoint variables. Similarly to obtain the co-state variables by using Pontryagin's Maximum Principle stated in literatures [51,52,55], with the existence result the following theorem is stated: be the optimal control and S * ; C * q ; H * p ; C * v ; C * i ; H * u ; H * a ; M * u ; M * a ; R * ; H * t be the associated unique optimal solutions of the optimal control problem (17) with initial condition (18) and objective functional (19) with fixed final time T f (20). Then there exists adjoint function l * i � ð Þ; i ¼ 1; : : :; 11 satisfying the

Numerical results
In this section we have presented the numerical result we have obtained using the parameters value collected in Table 3 below. We have collected data from a variety of sources, and have compiled the values in the table for the convenience of the constructed model numerical simulations and to verify the analytical results.

Numerical simulations and discussions of the deterministic model (3)
In this section, a numerical simulation of the entire HIV/AIDS and COVID-19 co-infection model given in Eq (3) is performed. We used ode45 fourth order Runge-Kutta scheme to examine the effect of various parameters on the spread and control of COVID-19 mono-infection, HIV/AIDS mono-infection, and HIV/AIDS and COVID-19 co-infection. The parameter values presented in Table 3 are used for numerical simulation. Moreover, we have investigated the stability of the endemic equilibrium point of the co-infection model (3), the effects of parameter on reproduction numbers, and the impact of treatment primarily on dually-infected individuals in the community.  Table 3. We can deduce from the figure that after a year, the solutions of the COVID-19 and HIV/AIDS co-infection dynamical system (3) are approaching the endemic equilibrium point of the given dynamical system whenever the co-infection effective reproduction number

Numerical simulation to show the effect of k 3 on R HM
The effect of the HIV protection rate on the HIV/AIDS effective reproduction number R HM is depicted in Fig 4. The graph shows that as the value of protection rate k 3 increases, the effective reproduction number R HM decreases and for k 3 > 0.771 indicates that R HM is reduced to less than one. As a result, the public health and policymakers must focus on increasing the values of the HIV/AIDS protection rate k 3 in order to control HIV/AIDS spread which may causes for existence of co-infection in the community.

Simulation to show the effect of κ on R C
A numerical simulation in order to show the effect of COVID-19 treatment on the COVID-19 effective reproduction number R C is given by Fig 5. The graph shows that as the value of the treatment rate raises, the COVID-19 basic reproduction number decreases and for the value of κ > 0.776 implies that R C < 1.  The graph shows that when the value of k 4 grows, the COVID-19 effective reproduction number decreases, and values of k 4 > 0.9 suggest that R C < 1: As a result, public health authorities must focus on increasing the COVID-19 immunization rate k 4 in order to prevent and control COVID-19 spread in the community. Biologically, this indicates that the COVID-19 infection reduces as the immunization rate k 4 rises.

Numerical simulations of optimal control strategies
To verify the analytical results, the optimal control model system (17) is simulated using the parameter values given in Table 3 with positive weight constants w 1 = w 2 = w 3 = w 4 = 18. The optimal control system is composed of two dynamical systems, the state dynamical system (17) and the adjoint dynamical system (27), each with its own initial and final-time conditions, with the control value state in Eq (26). The fourth forward-backward Runge-Kutta iterative method is used to solve this optimality system. The state Eq (17) is solved with the initial values of state variables using the fourth-order forward Runge-Kutta method. We used backward fourth order Runge-Kutta to solve the adjoint equations once we had the solution of the state functions and the value of optimal controls. To determine the impact of control measures on the reduction of the HIV/AIDS and COVID-19 co-infection we have the following three cases of optimal control strategies:

5.19.
Co-infection simulation with strategy 1 (u 1 ¼ 0, u 2 6 ¼ 0, u 3 6 ¼ 0, and u 4 6 ¼ 0) In this subsection simulation is done for the cumulated HIV/AIDS and COVID-19 co-infection when there is no control strategy in place and when there are controls involving COVID-

5.21.
Co-infection simulation with strategy 3 (u 1 6 ¼ 0, u 2 6 ¼ 0, u 3 ¼ 0, and u 4 6 ¼ 0) In this subsection simulation is done when there is no control strategy in place and when there are controls involving HIV protection, COVID-19 protection, and HIV treatment without COVID-19 treatment measure. Fig 21 shows the result that all the prevention and control strategies except HIV treatment strategy efforts are implemented, the number of individuals co-infected with HIV and COVID-19 decreases drastically to zero after 7 years.

5.22.
Co-infection simulation with 4 (u 1 6 ¼ 0, u 2 6 ¼ 0, u 3 6 ¼ 0, and u 4 ¼ 0) In this subsection simulation is done when there is no control strategy in place and when there are controls involving HIV protection, COVID-19 protection, and COVID-19 treatment without HIV treatment measures. Fig 22 shows the result that all the prevention and control strategies except HIV treatment strategy efforts are implemented, the number of individuals coinfected with HIV and COVID-19 decreases drastically to zero after 8 years.

5.24.
Co-infection simulation with strategy 6 (u 1 6 ¼ 0, u 2 6 ¼ 0, u 3 ¼ 0, and u 4 ¼ 0) In this subsection simulation is done when there is no control strategy in place and when there are control strategies involving protection strategies for COVID-19 and HIV single infection without HIV and COVID-19 treatment measures. Fig 24 shows the result that protection strategies efforts are implemented without treatment strategies, the number of individuals coinfected with HIV and COVID-19 decreases drastically to zero after 8 years later.

5.25.
Co-infection simulation with strategy 7 (u 1 6 ¼ 0, u 2 6 ¼ 0, u 3 6 ¼ 0, and u 4 6 ¼ 0) In this subsection simulation is done when there is no control strategy in place and when there are all the control strategies involving protection and treatment for both COVID-19 and HIV single infections. Fig 25 shows the result that all the protection and treatment strategies efforts are implemented, the number of individuals co-infected with HIV and COVID-19 decreases drastically to zero after 3.

Conclusions
In this paper, we formulated and investigated a continuous time dynamical model for the transmission of HIV/AIDS and COVID-19 co-infection with protection and treatment strategies. The mode incorporate four non-infectious groups the susceptible group, the HIV

PLOS ONE
protection group, the COVID-19 protection group, and the COVID-19 vaccinated group and this made the model highly nonlinear and challenging for the qualitative analysis of the coinfection model. The model has been mathematically analyzed both for the sub-models associating the cases that each disease type is isolated and in the case when there is co-infection. In addition an optimal control problem model that minimizes the cost of the infection as well as minimizes the control efforts to control the diseases transmission in the community is formulated and analyzed. The model includes the intervention strategies, protective as well as treatment and numerical simulations of both the deterministic model and optimal control problem models are presented. In the analysis it has been indicated that the effect of protection as well as treating the infected ones with the available treatment mechanisms affects significantly the optimal control strategy and its outcome. From the optimal control problem simulation results it can be concluded that applying both protective and treatment control mechanisms at the population level yields both economic as well as epidemiologic gains. Therefore, we recommended to the stake holders to give more attention and the overall optimal effort to implement both the protective as well as treatment control strategies to minimize the single infections as well as the co-infection diseases transmission in the community. This study did not considered the stochastic method, fractional order method, impacts of the environment, structure of human age, the spatial structure, and real population primary epidemiological data. Based on these limitations potential researcher can consider to extend this study.